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Abstract 

Despite the discovery of a Higgs boson h{ 125 GeV) at the LHC Run-1, its self-interaction has fully 
evaded direct experimental probe so far. Such self-interaction is vital for electroweak symmetry 
breaking, vacuum stability, electroweak phase transition, and Higgs inflation. It is a most likely 
place to encode new physics beyond the standard model. We parametrize such new physics by 
model-independent dinrension-6 effective operators, and study their tests via Higgs pair production 
at hadron colliders. We analyze three major di-Higgs production channels at parton level, and 
compare the parameter-dependence of total cross sections and kinematic distributions at the 
LHC (14TeV) and pp(100TeV) hadron collider. We further perform full simulations for the di- 
Higgs production channel gg —l hh —1 6677 and its backgrounds at the pp (lOOTeV) hadron 
collider. We construct four kinds of benchmark points, and study the sensitivities to probing 
different regions of the parameter space of cubic Higgs interactions. We find that for one-parameter 
analysis and with a Sab” * 1 (30ab _1 ) integrated luminosity, the gg —» hh —► 6677 channel can 
measure the SM cubic Higgs coupling and the derivative cubic Higgs coupling to an accuracy of 
about 13% (4.2%) and 5% (1.6%), respectively. 
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1. Introduction 


The LHC discovery of the light Higgs boson /i (125 GeV) [Tj has become a historical turning point 
of particle physics. The standard model (SM) |2j could provide such a Higgs boson [3], which joins 
three types of fundamental interactions: (i) the gauge interactions mediated by spin-1 weak gauge 
bosons ( W ., Z ); (ii) the Yukawa interactions with fermions mediated by the spin-0 Higgs boson h ; (iii) 
and the cubic and quartic Higgs self-interactions h 3 and h 4 . But the type-(ii) and type-(iii) Higgs 
interactions are largely untested so far, which provide the most likely place to encode new physics 
beyond the SM. The current ATLAS and CMS measurements [4] find the Higgs boson /i (125 GeV) 
to appear SM-like, but only have weak sensitivities to hrf and hbb Yukawa couplings, while even 
the LHC run-2 could not sensitively probe most of other Yukawa couplings via direct detection [5|. 
Furthermore, the LHC has little sensitivity to probing the type-(iii) Higgs self-interactions. It was 
shown that the high luminosity LHC(14TeV) with an integrated luminosity of Sab" 1 could probe 
the h 3 coupling to about 50% accuracy [6j 7], an d the improved analysis could reach a sensitivity 
about 30% — 20% [§]. With the current measurements of Higgs and top quark masses, the SM 
Higgs vacuum becomes unstable around 10 9 " n GeV [9] and is very sensitive to new physics [ TO] . 
So new physics is expected to enter the Higgs potential and modify its self-interactions well below 
Planck scale mm- The Higgs self-interactions are vital for the spontaneous electroweak symmetry 
breaking [2], the electroweak phase transition HD, and the Higgs inflation [12]. Hence, it is important 
to probe Higgs self-interactions with precision and pin down the associated new physics deviations 
from the SM. 

The SM Higgs sector is described by the gauge-invariant renormalizable Higgs potential, 

V = -n 2 H ] H + \(H ] H) 2 , (1.1) 


where H = (n + , ^(u + h + is the Higgs doublet and v — 246GeV denotes the vacuum 

expectation value (VEV). Thus, the Higgs self-interactions take the form, 


hint 


^ 3 + 


4k h 4 

4! n ' 


( 1 . 2 ) 


where at tree-level we have the cubic and quartic couplings of the SM Higgs boson, A 3 = 6Au = 
‘IMl/v and A 4 = 6A = 3 M^/v 2 . Hence, given the observed Higgs mass Mh ~ 125 GeV [1], the Higgs 
self-couplings are completely determined in the SM. One could naively make a shift of Higgs coupling 


within the SM Higgs potential (1.1), A —> A 7 = A + 5 \, but it causes no observable effect, because 
this just redefines the renormalizable Higgs coupling as A 7 no matter what value 5 A would take. 
The nontrivial modification of Higgs couplings could only arise from higher dimensional effective 
operators whose effects cannot be absorbed into the dimension-4 SM Lagrangian. 
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Given that the SM Lagrangian already contains all possible gauge-invariant and renormalizable 
operators up to dinrension-4 and no new physics is found yet, the possible leading new physics devia¬ 
tions from the SM can be generally parametrized by gauge-invariant dimension-6 effective operators 
in a model-independent way m- Since the Higgs potential acts as the core of spontaneous elec- 
troweak symmetry breaking and has escaped from direct measurement so far, it stands out as a most 
likely place to encode new physics beyond the SM. Such new physics will certainly modify the Higgs 


self-interactions (1.2), via dimension -6 operators, which may not only shift the Higgs self-coupling 
itself [due to the operator (H^H) 3 ], but also modify the structure of Higgs self-interactions (due 
to the dimension -6 derivative operators). To sensitively probe such new physics in the cubic Higgs 
self-coupling, it is important to study di-Higgs production at high energy hardon colliders [331 [15] Q 

For hadron colliders, the main di-Higgs production channels include gluon fusion production, 
vector boson fusion (VBF) production, and top-pair associated production. Over a wide energy 
range, the total cross section of di-Higgs production via gluon fusions is almost 10 times larger than 
the other channels mm- Hence, it provides the dominant di-Higgs production. The decay mode 
hh —> 6677 has much cleaner background than others, so it has attracted efforts from both theoretical 
and experimental sides [ 6 ] [18] [19j for studying the potential of the high luminosity LHC (14TeV) and 
the future pp (lOOTeV) collider. Other di-Higgs decay modes with larger signal rates are also explored, 
such as hh — > bbrr , hh — > bbWW* — > bb2£2v , and hh —»• bhbh, etc [20]. Due to large backgrounds in 
these channels, more elaborated strategies like boosted kinematics are needed. Another decay mode 
hh —> WW*WW* —> 3i3vjj was considered with the use of m T2 observable [21] . Some rare final 
states were also explored for pp(100TeV) collider [22]. In addition, two more production channels 
have received recent attentions. The top-pair associated production pp —>• tthh turns out to be 
complementary to gluon fusion gg — > hh with hh —?- bhbh final states [23]. The VBF production 
channel pp — > hhjj receives a large contribution from gluon fusion production in the signal region, 
which makes the VBF contribution almost negligible [23] - Most previous studies for the di-Higgs 
production focused on the SM Higgs potential. There are recent analyses studying the contributions 
of dimension -6 operators to gg —>■ hh with hh —> 6677 [25j and hh —> bbrr [2B], It was noted that 
certain new operators can modify kinematic distributions of the final states as well as total cross 
section. For an operator that induces tthh coupling, the kinematics could be useful to increase the 
sensitivity [27, 28j- In general, including these operators with associated new coefficients will enlarge 
the parameter space of new physics, and thus make the probe of each individual parameter in the 
cubic Higgs interaction harder. Certain simplifications are needed to reduce the large parameter 
space. 

In this work, we will systematically analyze the new physics contributions of dimension -6 operators 
to the di-Higgs productions. For good physics reasons, we will focus on two rather unique bosonic 
dimension -6 operators which contribute to the cubic Higgs coupling and build a 2-dimensional (2d) 


1 Measuring Higgs quartic coupling would be even more challenging in the foreseeable future m- 
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parameter space. In particular, we will inspect the new operator that induces derivative cubic Higgs 
coupling, and thus has enhanced contributions to high energy processes. We will derive nontrivial 
perturbative unitarity constraints on these dinrension -6 operators. Then, we study the di-Higgs 
production via three major channels for probing cubic Higgs couplings. For this, we will perform a 
parton level analysis at the LHC (14TeV) and pp (lOOTeV) collider. Finally, we present a full analysis 
(including Delphes 3 fast detector simulations) for the di-Higgs production gg —> hh with hh —> 6677 
at the pp (lOOTeV) collider. From this we study the probe of new physics scales associated with the 
dimension -6 operators. We also find nontrivial interference between different operators, which can 
be probed by using relevant kinematic distributions. 

This paper is organized as follows. In section[2j we discuss the dimension -6 operators relevant to 


Higgs self-interactions, and identify the unique operators (2.16) which spans a 2d parameter space. 
We also motivate these operators by nonminimal Higgs-gravity interaction. We further study the 
perturbative unitarity constraints on the cutoff scales associated with dinrension -6 operators. In 
section]!] we analyze three major di-Higgs production channels at parton level and compare the 
parameter-dependence of total cross sections and kinematic distributions. In section]!] we perform 
full simulations for gg —>• hh —>• 6677 at the lOOTeV hadron collider, and study the sensitivity to the 
2d parameter space for four benchmarks. We conclude in section]!] Finally, Appendix A discusses 
the redundancy of dimension -6 operators, and Appendix B summarizes the loop functions of triangle 
and box diagrams for the analyses of sections[3]j4] 


2. New Higgs Self-Interactions from Dimension-6 Operators 

2.1. Identifying Relevant Dimension-6 Operators 


The SM Lagrangian is a fairly good effective theory up to gauge-invariant renormalizable operators of 
dimension-4. The possible leading new physics deviations are generally parametrized via dimension -6 
effective operators]^] 


£ 


eff — 



(2.3) 


where A characterizes the cutoff scale, and the dimensionless coupling f n is expected to be around 
0(0.1 — 1) for each given operator (unless suppressed by extra symmetry). The LHC Run-1 data [3j 
have constrained the 125 GeV Higgs boson to be fairly SM-like and found no new light particle beyond 
the SM. Hence, it is well-motivated to use the standard effective theory formulation of possible new 
physics effects via dimension -6 operators m, and assume that no other light field exists below its 

2 If the light neutrinos turn out to be Majorana fermions, there is a unique dimension-5 effective operator 1291 . 

which provides Majorana neutrino masses and violates lepton number by two 
units. Here C = * 7 2 7 ° is the charge-conjugation operator, (i,j) are flavor indices of left-handed lepton doublet, 
and (a, a', /3, /3') are indices of SU(2) doublets. This dimension-5 operator is irrelevant to our current study of the 
Higgs self-interactions. 
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cutoff scale. The full set of gauge-invariant dimension -6 operators that modify Higgs self-interactions 
includes [30], 




o 


< 3>,3 




0*4 = {Di*H?(D IJi H)(HiH). 


(2.4) 


Among all four operators, 0$, 2 and 3 modify scalar sectors only, while 1 and 4 also 
contribute to gauge boson masses and couplings. The operator O $ ] contributes to the mass m z , but 
not to m w . Thus, it violates the custodial symmetry and is severely constrained by the electroweak 
precision parameter T. For collider searches, it is safe to neglect the effects of 0$ 1 [31]. With 
the equation of motion (EOM), there is redundancy among dimension -6 operators. As explained in 
AppendixjAj the subset operators (0$ 2 , 0$ 3 ) 4 ) hi (2.4) are not independent. Including the SM 

Yukawa interactions, another type of dimension -6 operators O $ j become relevant, 


0 9>f = {H*H)LHf R + h.c., 


(2.5) 


where L = (/£, f£) T denotes the SU(2)l doublet, and f R the SU(2) L singlet. Among all operators 
mentioned above, one operator can be eliminated via EOM. We choose to drop 4 hereafter. Thus, 
we have two rather unique bosonic dimension -6 operators (O# 2 , 0 $ 3 ) relevant to the present study 
of Higgs self-couplings. 

Next, we inspect the contributions of (£)$ 2 , O $ 3 , O $ j) to the Higgs self-couplings, as well as the 
Higgs-gauge and Higgs-fermion couplings. For later convenience, we define a dimensionless coefficient 


Xj and an effective cutoff scale A ; - for each operator in (|2.4|)-(|2.5|), 


/<t 


x i = , 


Aj - 


A 


l/< 


( 2 . 6 ) 


EH 


Since (D§ 3 is a non-derivative operator, it only affects Higgs mass and self-couplings. In par¬ 
ticular, it modifies the relation between the observed Higgs mass and cubic Higgs coupling. The 
derivative operator 2 induces the following term for Higgs field, 


O 


$,2 -t -^{h + vfd^hd^h, 


with x 2 defined in (2.6). This modifies the Higgs kinetic term as, 


£ ki n = \ (1 + * 2 ) d^hdph . 

Thus, we can dehne the canonical Higgs field via rescaling h —> (h , with the factor, 

C = (l + X 2 )-5. 


(2.7) 


( 2 . 8 ) 


(2.9) 


This induces a universal modification to all Higgs couplings with SM particles. After the normal¬ 


ization, Eq. (2.7) also generates a derivative cubic Higgs interaction, ^x^hd^hdji. In contrast to 
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the SM cubic Higgs coupling, this new derivative interaction vertex will be enhanced by the center of 
mass energy in high energy processes, and thus may have distinctive kinematic feature. The modified 
cubic Higgs coupling is 


h-h-h: -iM C ,i- 

v 

.c 


X3(2 3^i) 


= -i— [3 (1 + f) Ml - x (pj +p\+ pi) ] 


( 2 . 10 ) 


In the above, M h is the physical mass of the Higgs boson, which receives contributions from both 


the kinetic rescaling factor (2.9) and the dimension-6 operator O $ 3 . So, we deduce the Higgs mass 
formula, M% = M| 0 [ 1 — x 3 /(2A)]£ 2 , where M h0 = y/2Xv is the SM Higgs mass. We see that M h 
depends on {£, x 3 }. For convenience, we replace (x 2 , x 3 ) by another two independent inputs (f, x) 
which parametrize the modifications of cubic Higgs coupling with different kinematic properties, 

2v 2 


r = -x 3 C 


3 M 2 ’ 


x = x 2 ( . 


( 2 . 11 ) 


With x, the rescaling factor can be rewritten as ( = (1 — x) 1 / 2 . We also note that the operators 
2 , 0$ 3 ) do not affect the W mass at tree-level, so the Higgs VEV is determined by the Fermi 
constant G F as in the SM, v = ( V2G F ) ~ 246 GeV. The modification to Higgs-gauge boson 

coupling only arises from rescaling the Higgs field, 


Vp-Vu-h: 


V^-V u ~h-h: 


v 


(2.12a) 

(2.12b) 


where V = W, Z. In unitary gauge, the Higgs-fermion dimension-6 operator (2.5) generates following 
term, 


° tj ^/k^ (v+h)3ff - <2 ' 13) 

This contributes to fermion mass, rrij = 

At the same time, it modifies the / — / — h Yukawa coupling. Replacing yj 11 by rrij, we deduce the 
following effective Yukawa coupling, 


V2 


_ 1 x ^j ( where yj 11 is the SM Yukawa coupling. 


f-f-h: 


— 1 - 


. Cm, 


1 - 


x f v 

y/2 m. 


This operator also induces a dimension-5 vertex h 2 ff , 

- 0 3X r 

f-f-h-h: iC 2 f 


V2v ' 


(2.14) 


(2.15) 


It contributes to the gluon fusion production gg —> hh with triangle quark loop. Since top quark is 
most relevant in practice, it is natural to set f = t for the present analysis. 
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The dimension-6 operators (2.4) and (2.5) are subject to constraints from measurements of single 


Higgs production at the LHC. The current data put the best bound on Higgs-gauge couplings (2.12a) 
m- For a future e + e Higgs factory with 250 GeV collision energy, the sensitivity to e + e —> Zh 
cross section is expected to be 5a/a = 0(0.5%) [32] with a 5ab _1 integrated luminosity. This is 


a direct probe of the modification of Higgs-gauge couplings and thus constraints |x| at 1% level 
J. The operator 0$ 3 will contribute to the e + e - —> Zh cross section via one-loop corrections 
From this, the sensitivity to A| m is estimated to be about 35% at the e + e~ Higgs factory with 
a 5ab _1 integrated luminosity. Besides, many other dimension-6 operators can contribute to the 
gauge boson kinetic terms and thus the wavefunction renormalization. This will further shift Higgs- 
gauge couplings, and make the constraint on individual operators much weakened. For top Yukawa 


coupling, the LHC run-2 has weak sensitivity to probing the deviations in (2.14). The precision of a 
high-luminosity LHC (HL-LHC) is expected to be around 10% 


Dihiggs production at high energy hadron colliders is an important way to measure the cubic 


Higgs coupling. The dimension-6 operators (2.4)-(2.5) contribute in different di-Higgs production 
channels. For gluon fusion and top-pair associated production, three operators 0$ 3 , 0$ 2 and t 
are relevant. For the two operators that modify Higgs self-interactions, 3 contributes to the SM 
cubic Higgs coupling by a simple shift (without affecting its Lorentz structure), and is commonly 
studied in the di-Higgs production literature. On the other hand, O $ 2 induces derivative cubic 


Higgs coupling (2.10), and is rarely studied for di-Higgs production. This operator contributes to 
the Higgs-gauge coupling via Higgs wavefunction renormalization and thus may receive constraint 
from measuring single Higgs productions (via vector boson fusion or Higgs-gauge-boson associated 
production) at colliders. But, since other dimension-6 operators also contribute to the Higgs-gauge 
couplings with interferences and possible cancellations, there is no unique constraint on O ^ 2 J^] 
Hence, it is important to directly probe the derivative cubic Higgs coupling induced by O $ 2 via 
di-Higgs production, which has distinctive kinematic features from other non-derivative operators. 
For the present study, we will focus on the new physics contributions to the Higgs self-couplings in 
di-Higgs production, and drop the fermionic operator O^ t (which was considered before [2?] [2S| and 
is irrelevant to Higgs self-interactions Q With these considerations, we define our parameter space 
by identifying the two rather unique bosonic dimension-6 operators, 


0*, 2 = \d^(H^H)d^H), 0* t3 = \(H^Hf. (2.16) 

3 One could expect other possible precision constraints on 0$ 2 from such as the muon anomalous magnetic moment 
g^ — 2 at two-loop. Again, other new physics operators such as the dimension-5 Pauli term Fj t iY cr , “''0 (with if> being 
muon field) can also contribute to g^— 2 at tree-level and become dominant. Hence, there is no unique constraint on 
2 from g^—2 at two-loop. 

4 In principle, 0 ®, t could be discriminated from 0® 2 and Cl ®,3 by further performing a combined analysis of three 
di-Higgs production channels via gluon fusion, VBF production, and top-pair associated production. With these and 
the single Higgs production gg —> h , we may also discriminate another operator G al ‘ 1 ' G'*H' H (which does not modify 
the Higgs self-coupling). It is possible that some other dimension -6 operators may contribute to the backgrounds as 
well, but without any special cut or selection they are expected to be much smaller than the SM backgrounds (from 
the dimension-4 operators of the SM). For clarity of the current analysis, we assume that these additional operators 
are negligible. 
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In the following subsection 2.2 we will further motivate the operator 2 from the Higgs-gravity 


interaction. Then, we derive generic perturbative unitarity bound on 2 in Sec. 


2.3 


2.2. Motivation from Higgs Gravitational Interaction 


The world is apparently described by a joint effective theory of the SM and general relativity (GR) 
up to accessible energy scales so far. It is important to probe the interface between the SM and GR. 
With the LHC discovery of a light Higgs boson h (125GeV), there is a unique dimension-4 operator 
at this intersection, namely, the nonminimal interaction between the Higgs doublet H and the Ricci 
scalar-curvature 77 


S t = 


jd 4 x y/—g H1Z , 


(2.17) 


where £ is a dimensionless coupling. With the proper normalization of graviton propagator, it is 
clear that under perturbative expansion the coupling £ is always associated with the suppression 
factor 1/Mpj. Hence, £ 1 can be well consistent with perturbative calculation. The current 

LHC constraint on this coupling is actually rather weak, and £ can be as large as O(10 15 ) [37) [38]. 
Nontrivial constraints from perturbative unitarity were derived before f3S|. The operator (2.17) 
has many physical applications such as the Higgs inflation m, gravitational dark matter [39] . and 
collider signatures [38J. Including this operator, we write the joint effective Lagrangian of the SM 


and GR, 


Sj = 




-g( J ) 


\m 2 + £tft h] 77^ - J2l F ^ F r + (D^iD^H) - V{H) 


(2.18) 


where 77*^ is the Ricci scalar corresponding to the Jordan frame metric g\?J, and Fj* wi = (W“„, B^ v ) 


are gauge field strengths of the electroweak gauge group SU(2) L <8> U{1) Y - In (2.18), we can readily 
include the SM fermionic Lagranian C F as well, though it is not relevant to the discussion below. For 
practical applications, it is convenient to make a Weyl transformation for metric field, ()\\y = , 

with the factor 


„ M 2 + 2fH^H 

fr = - 


M p\ 


(2.19) 


After changing variable, we write down the action with new metric g)w\ 


S E — 


I 


- E + sips (S^H)) 


1 Mln - . 

2 p ' 4 

3 




( 2 . 20 ) 


For simplicity, we drop the superscript (E) for all geometric quantities associated with g FU here. 
Since the nonminimal interaction term is transformed away and the gravity sector becomes normal, 
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the new metric is called Einstein frame. In this case, all effects of £ appear in matter sector and are 
represented by a series of higher dimensional effective operators. Expanding these ^-induced terms to 


leading order, we can deduce two relevant dinrension-6 Higgs operators 2 and 3 from (2.20), 

( 2 . 21 ) 


a §2 


associated with two different cutoff scales, 

Mpi 


v = 


£ 


— 


Mp\ 

vr' 


( 2 . 22 ) 


is 


Among dinrension-6 operators in (2.4), A^ is related to 9 with / $9 /A 2 = 6/A^, which i 
generated due to the third term of Eq. (2.20). Expanding the l/fl factors in (2.20) will induce 
(0 $ 3 , 4 , j), with a cutoff characterized by A^ 2 = Mpi/y^ . For the operator O $ 3 , we have 

/,j> 3 /A 2 = 12/A| 2 . The other two operators 0$ 4 and j are induced from 1/fl expansion with 
the following coefficients, 


A 2 


Off. 4 + 


A| 2 


a 


<£>,/ ’ 


(2.23) 


where is the SM Yukawa coupling of the fermion /. The effective theory with such Higgs-gravity 
interactions can be viable for a wide range of £. To be relevant to collider physics, we need £ S> 1 
[380 which implies A| 4 -C A| 2 . Hence, in this effective theory, the operator O ^ 2 will give dominant 
contributions, while other operators 3 and 4 , O ^ j) are negligible. 


2.3. Constraints from Perturbative Unitarity 


In this subsection, we derive perturbative unitarity bound on the parameter space of dimension- 


6 operators defined in (2.16). We analyze the longitudinal weak boson scattering and top-Higgs 
scattering in high energy regime. We find that their scattering amplitudes are largely enhanced by 
E 2 and E 1 contributions from the derivative cubic Higgs couplings, and would eventually violate 
perturbative unitarity with the increase of scattering energy. This places an upper bound on the 
validity range of perturbation expansion of the effective theory, above which certain nonperturbative 
dynamics or new physics have to set in0 For the current analysis, we will derive perturbative 
unitarity bounds for both types of processes. Since the energy dependence of gg —>• hh amplitude is 
rather mild, it cannot place better bounds than the processes mentioned above, and thus needs no 


consideration here. 

5 As we clarified before [38], in this effective theory formulation, we do not concern any detail of the UV completion 
above the cutoff 2 - There are many well-motivated TeV scale quantum gravity theories on the market. For 
instance, extra dimensional models with compactification scale at = O(lOTeV) will reveal the Kaluza-Klein modes 
at energies above this scale, and other related UV dynamics may show up above this cutoff as well. 

6 Since the joint effective theory of SM+GR is nonrenormalizable and its UV completion is unknown, any naive partial 
resummation within this effective theory itself cannot give reliable unitarity restoration [38]. Hence, the perturbative 
unitarity bound is important for such nonrenormalizable effective theories. 
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Figure 1: Longitudinal weak boson scattering processes, VlVl —> hh (VlVl), where V = W ± , Z°. 
The crossing channels also give gauge-Higgs boson scattering. 


Fig.[l] depicts the longitudinal weak boson scattering VlVl —> hh (VlVl) and the gauge-Higgs 
boson scattering in the crossing channels. The new physics of dimension-6 operators modifies the 
Higgs-gauge coupling and the Higgs self-couplings, which can induce nonzero 0(E 2 ) enhancement 
in the scattering amplitudes [38]. Since dimension-6 operators are gauge-invariant, the longitudinal- 
Goldstone boson equivalence theorem (ET) [3D] can be established [38] . Hence, the same E 2 en¬ 
hancement must show up in the corresponding Goldstone boson scattering amplitudes. To derive 
the optimal unitarity constraints on dimension-6 operators, we perform a coupled channel analysis 
of all electrically neutral channels for Goldstone boson and Higgs boson scatterings, with initial/final 
states {|7r + 7r - ) , |7r°7r°) , 4^ | hh) , |7T°/i)}. We compute the relevant leading scattering amplitudes 

at 0(E 2 ), 


T[VT + 7T _ —> 7T + 7T~] 
T[7T + 7T~ —> 7T 0 7r°] 

T[7r + 7r _ —> hh] 
T[ir°h -> 7 T°h] 

E[tT°TT°—> 7T 0 7r°] 


(l+cos0)E 2 


2v 2 


E 2 


T[7T U 7T U — > hh] = x(\ — x) 


E 2 


—X (1 — 2 ?) 


(1 —cos 0)E 2 
2V 2 


O(E 0 ), T[hh ->• hh] = O(E 0 ), 


(2.24) 


where E is the center-of-mass energy and 8 denotes the scattering angle. With these, we compute 
the corresponding partial wave amplitudes, 


a e (E) 


1 f 1 

- / dcos 8 Pp(cos 6)T(E, 8). 

327r 


(2.25) 


We perform a coupled channel analysis for the in/out states {|7r + 7r ), ^ |7r°7r°) , ^ | hh) , |7r°/i)}, 
Then, we can derive the following 4x4 matrix for the s-wave amplitudes at 0(E 2 ), 


a o(E) 


< 1 

xE 2 \/2 

327 TV 2 y/2(l-x) 

V 0 


\J~2 V2(l-x) 

0 1-x 

1 — x 0 

0 0 


0 N 

0 

0 

-(1-x) j 


(2.26) 


For a sizable |1 — x|, the scattering amplitudes with Higgs in initial/final states have dominant 
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contributions. We deduce the following eigenvalues, 

xE 2 


a, 


o 


dms (E) = 32^2 diag(^l + v / l+3(l-x) 2 , 1- ^\+2,{l-x) 2 , -(1-x), -l) , 


(2.27) 


and impose the s-wave unitarity condition |Rea 0 | < ^ on the maximal eigenvalue. Thus, we derive 
the perturbative unitarity bound on the scattering energy, 

\/ l 67 tv 


E < A U;L — 


(2.28) 


[|x|(l + Vl+3(l-x) 2 )] 1 /2 ‘ 

We plot this bound A u;l as a function of x in Fig.[2](a) , where the blue region (including the overlap 
with red region) denotes perturbative unitarity violation. We also show the dependence of unitarity 
bound on the effective cutoff A 2 of the dimension -6 operator 0$ 2 in plots (b) and (c) for x 2 > 0 
and x 2 < 0, respectively. For small |x|, we find A U1 ~ ^16^/3 A 2 at leading order. As mentioned 
earlier, x could be constrained by measurements of Higgs-gauge coupling in single Higgs production 
due to its contribution to the rescaling of Higgs kinetic term. But, given the contributions from 
other dimension -6 operators to the Higgs-gauge coupling and their possible large cancellations with 
that of 2 , the Higgs-gauge coupling could be SM-like while x 2 is more or less free from this 
constraint. In this case, O# 2 still receives general perturbative unitarity bound from high energy 
scattering processes involving its induced derivative Higgs self-couplings, even though Higgs rescaling 
effect may be negligible. Thus, we derive the corresponding unitarity bound by turning off the Higgs 


rescaling effect in (2.26), 


E < A' m = 


\J 16-7T \ 


(2.29) 


s 1 / 4 !^ 1 / 2 ' 

We depict the upper bound (2.29) by the blue dashed curve in Fig.j2|a)-(c). We see that A(jj turns 
out to be weaker than the bound Api- In the later analysis of di-Higgs production via vector boson 
fusion, we will be conservative and select signal events by imposing the weaker bound < A^j. 

Fig. I presents Feynman diagrams for tt —> hh(V L V L ) scattering, where V = W ±, Z°. In high 
energy limit, the leading amplitudes from dimension -6 operator 0^ 2 are enhanced by E 1 terms. 
According to equivalence theorem, we compute the leading amplitudes with final state V L V L replaced 
by the corresponding Goldstone bosons. Among all contributions, the amplitudes with i/it-channel 
quark-exchange and the SM Yukawa coupling approach constant in high energy limit. Only the s- 


channel Higgs-exchange with cubic derivative Higgs coupling in (2.10) gives the O^E 1 ) asymptotical 


behavior and may violate perturbative unitarity. To derive the optimal bound, we define the spin-0 

1 Y 

and color-singlet helicity state of top-quark pair, i.e., \tt) s = > ) — |f“ )) 

v 2-^c Q=1 

we compute the scattering amplitudes at the leading 0(E 1 ), 


Thus, 


T[\tt ) s -> |tt + tt )] = T[\tt) 
T[\it) 


I 0 0\1 a2 m tE 

|vr u 7T u )J = — v 6 x C — Y~ 


| hh)\ = -\/ 6 xC 4 mt F 


T[\it) a -+ \7T°h)} = O(E 0 ) 


(2.30) 
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Figure 2: Perturbative unitarity violation region from weak boson scattering (blue) and top-Higgs 
scattering (red) as a function of x in plot-(a), and a function of A 2 in plots (b) and (c) for X 2 > 0 
and X 2 < 0 , respectively. 



Figure 3: Feynman diagrams for tt —> hh {V L V L ) scattering, where V = W±, 


Z°. 


where E is center of mass energy. To optimize the unitarity bound, we can further define an 0(4) 
singlet final state \S) = (2 |7r + 7r _ ) + 17r°7r 0 ^ + | hh)') . Hence, we derive, 


T[\tt) s -^ |S)] = —x(l — x)(4 — x) 


\/3m t E 

2 ^ 


(2.31) 


Using (2.25), we compute the partial wave amplitude and impose the s-wave unitarity condition 
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|Rea 0 | < ^ . With these we deduce the perturbative unitarity bound on scattering energy E , 


E < A U2 — 


167T7r 


\/3 m f |x(l-x)(4-x)| 


(2.32) 


We plot the upper bound A U2 in Fig. [2] (red contours) as a function of x and A 2 , respectively. For 
small |x|, we derive A U2 ~ A-kA^/ V3m t at leading order. It is clear that the bound from top-Higgs 
scattering is much weaker than that of the weak boson scattering. 


3. Higgs Pair Production at Hadron Colliders 


In this section, we study di-Higgs production for the effective theory defined in Eq. (2.16) at both the 


LHC(14TeV) and future pp(lOOTeV) collider. There are two new parameters (x 2 , x 3 ), which may 


be reparametrized as (x, f) in Eq. (2.11) for convenience. The major di-Higgs production channels at 
high energy hadron collider include gluon fusion production ( gg —> hh), top-pair associated produc¬ 
tion (pp —> tthh ), and VBF production (pp —> hhjj). In the following, we analyze these production 
channels at parton level, and compare their differences in total cross sections and in kinematical 
distributions over the parameter space of (x, r). 


With the modified cubic Higgs couplings (2.10) from dimension-6 operators, we derive the differ¬ 


ential cross section for gluon fusion production, 


dd(gg -t hh) 
dt 


G 2 F a 2 s 

512(2vr) : 


c 4 


(1 +r\ 


3m,i 


s — mt 


— x 


s + 2nij l 


>h s ~ m h 

where (s, t) are partonic Mandelstam variables, and (F^ G n ) are loop functions given in Ap- 
pendix[Bj The new contributions from x 2 (x) arise in two ways. The first is an overall rescaling 
factor C 4 of the cross section, and the second is contributed by the derivative cubic Higgs coupling. 
The parameter x 3 only appears in r , which shifts the SM cubic Higgs coupling. We generate signal 
events by MadGraph5 pEjJ^] The QCD corrections can be significant [43], but they are insensitive 
to the structure of cubic Higgs coupling]^] so we normalize the cross section at (f, x) = (0, 0) to the 
SM NLO prediction m and implement the same A'-factor for full parameter space of (r, x). For 
gluon fusion, we have K = (2.27, 1.44) for \/s = (14, 100) TeV. But, for analyzing the ratio of the 
cross section over that of the SM, it is rather insensitive to the A'-factor. We perform numerical fits 
for the total cross sections over the range — 1 ^ r ^ 1 and — 1 ^ x ^ 0.5 at both LHC (14TeV) and 
pp(100TeV) collider, 


F a + F c 


+ IQdI" 


(3.33) 


<7(59 - 

a hh) 

°(gg - 

hh) sm 


14TeV 


= (1-x) 2 (l - 0.83f + 3.7x + 0.29r 2 + 4.2x 2 - 2.0fx) , 


°{gg 

A hh) 

°(gg - 

hh) sm 


100TeV 


= (1-x) 2 (1 - 0.72 r + 3.6x + 0.22r 2 + 4.3x 2 - 1.7rx) . 


(3.34a) 

(3.34b) 


'To include the effect of finite top mass, we use the model file SMEFT_FF_bt for events generation. The relevant 
code is available at https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/HiggsPairProduction. 

8 As shown in Ref. [14] , for various dimension-6 operators relevant for gluon fusion production, the correction to the 
K -factor is around several per cent. 
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Figure 4: Cross sections of di-Higgs production via gluon fusion (blue), top-pair associated produc¬ 
tion (purple) and vector boson fusion (red) at the LHC (14TeV) (left plot) and pp(lOOTeV) collider 
(right plot). For each production channel, the (dashed, solid, dotted) curves depict cross sections as 
functions of x under three inputs of r = (—1, 0, 1). 


This shows that the fitted cross section ratio is not sensitive to the variation of collision energy 
from yfs = 14TeV to yfs = 100 TeV. This is mainly due to the m^/s suppression in the loop 
functions F A and under high energy limit [cf. Eq. ( |B.52| )]. Expanding (13.331) around the SM 
values (r, x ) = (0, 0), we derive the r dependence, d(a/a sm )/dr ~ —(0.7 —0.8). For the parameter 
x , the prefactor (1 — x) 2 = C 4 in Eq. (3.34) comes from rescaling factors of Higgs fields hh in 
the final state, while x in the second parentheses is contributed by the derivative cubic Higgs 
coupling. We note that these two contributions have some cancellation. For x > 0 (x < 0), 
the contribution from derivative coupling interferes constructively (destructively) with the SM part 
of (r, x) = (0, 0), while the total cross section is suppressed (enhanced) by the Higgs rescaling 
factor (1 — x) 2 . The blue curves in Fig. [4] depict the gluon fusion cross sections at pp( 14TeV) and 
pp(100TeV). The (dashed, solid, dotted) curves present the cross sections varying with x, under 
inputs r = (—1, 0, 1), respectively. From Fig.[4j we see that the di-Higgs production cross sections 
from gluon fusion exhibit a minimum in the x < 0 region, and the location of this minimum varies 
with the input value of r. 

In Fig.[5j using MadAnalysis-5 package [45], we present the normalized kinematic distribution of 
final state Higgs bosons at pp(100TeV) collider. The first column display the leading Higgs p T (h ) 
distributions; while the second column depict the M hh invariant-mass distributions of the Higgs pair. 
The shapes of distributions at the LHC(14TeV) and pp(100TeV) collider have some similarity since 
the cross section only has mild energy dependence. In the first row of Fig.[5j we have input r = 0, 
and the (blue, red, green) curves correspond to x = (—1, 0, 0.5); while the second row has x = —1, 
and (blue, red,green) curves correspond to r = (—1, 0, 1). For the parameter range x < 0, there is 
large cancellation between the SM box-loop diagram and the triangle-loop diagram with s-channel 
Higgs and new derivative cubic Higgs coupling over the intermediate momentum range. This makes 
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Figure 5: Parton level distributions of gg —> hh for the leading Higgs p T (1st column) and the 
invariant-mass M hh (2nd column) at pp(lOOTeV). In the first row, we input r = 0, and the 
(blue, red, green) curves correspond to x = (—1, 0, 0.5). In the second row, we input x = —1, 
and the (blue, red, green) curves correspond to r = (—1, 0, 1). 


the distribution more sensitive to r. In particular, if we turn off the SM cubic Higgs coupling by 
setting r = — 1, the events are mostly populated in large p T and M hh regions, as shown by the 
blue curves in the second row of Fig. [5} For x > 0 and r > — 1, all contributions add to each other 
constructively, and the normalized distributions do not significantly change^ 

Next, we consider the top-pair associated di-Higgs production. The dependence of its cross section 


on (x 2 , x 3 ) can be reparametrized in terms of (x, r), and is similar to that of (3.33). We generate 
the signal events by MadGraph5, and find the factor K = 1.2 for total cross sections at both the 
LHC(14TeV) and pp(100TeV) collider [ 171 . We perform numerical fits of total cross sections for 
—1 ^ r ^ 1 and —1 ^ x ^ 0.5, which are summarized as follows, 
a(pp —> tthh) 


a(pp - 
a(pp 


tthh) s 
■ tthh ) 


14TeV 


= (l-x) 2 (l + 0.23?- 0.73 x + 0.04r 2 + 0.60 x 2 — 0.26 rx), 


<j(pp —> tthh) s 


lOOTeV 


= (l-x) 2 (l + 0.23f- 0.80x + 0.07f 2 +2.2x 2 -0.54fx). 


(3.35a) 


(3.35b) 


9 In passing, Ref. [3S] studied interference between the SM cubic Higgs coupling and other SM contributions in a few 
di-Higgs production channels, with focus on the variations of collision energy and parton distribution function. 
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Figure 6: Parton level distribution of pp —> tthh for the leading p T distributions of Higgs boson 
(1st row), the invariant-mass distributions M hh (2nd row) at the LHC(14TeV) (1st column) and the 
pp(lOOTeV) (2nd column). In each plot, we set r = 0 , and input x = (—1, 0, 0.5) which correspond 
to (blue, red, green) curves, respectively. 

In comparison with the di-Higgs production via gluon fusion, the cross section of top-pair associated 
production is less sensitive to the change of either r or x , due to the dominance of diagrams irrelevant 
to Higgs self-interaction. But, the x-dependence of top-pair associated production cross section is 
much more sensitive to the increase of collision energy than that of gluon fusion production, especially 
for the x 2 term. We note that the derivative cubic Higgs coupling term interferes destructively 
(constructively) with the SM f/ti-channel exchange of top for x > 0 (x < 0). Hence, this process 
is complementary to gluon fusion production. In Fig. [4] we plot the total cross sections of top-pair 
associated di-Higgs production by purple curves. It is much suppressed in x > 0 region due to the 
overall rescaling factor (1 — x) 2 = £ 4 . For r > 0, it adds positive contributions to that of the SM, 
and makes the test of r easier [ 23] . 

We present in Fig. [6] the normalized kinematic distributions for top-pair associated di-Higgs pro¬ 
duction at parton level. The first row shows the leading p T distribution of the Higgs boson, and 
the second row depicts the di-Higgs invariant-mass (M hh ) distribution, at the LHC (14TeV) (in first 
column) and pp(100TeV) collider (in second column). At the LHC, they are rather insensitive to the 
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variation of (x, r). However, the pp(lOOTeV) collisions significantly improve the sensitivity to x. 
In comparison with the di-Higgs production via gluon fusion in Fig.[5j the top-pair associated pro¬ 
duction is more sensitive to the derivative cubic Higgs coupling, with more signal events populated 
in the higher p T and larger M, . region. To maintain perturbative unitarity, we will require signal 


events to obey Mhh < A U2 , where A U2 is derived in (2.32). We find that this bound at x = 0.5 is 
too weak to be relevant; and there are 77% (97%) signal events passed this requirement for x = — 1 
at y/s = 100 TeV (14TeV). 

Finally, we turn to the di-Higgs production via vector boson fusion, pp —> V*V*jj —» hhjj . 
Its cross section depends on (x 2 , x 3 ) through the overall rescaling factor £ 4 , the modified (SM-like) 
cubic Higgs coupling r , and the new derivative cubic Higgs couplings x . We generate signal events 
by Madgraph5 with electroweak process, and apply the following VBF cuts to two tagging jets m , 

14TeV: 2 < \rjj\ < 5, r/ ?| - ■q J2 < 0, p T - > 25GeV, M- > 500GeV; (3.36a) 

100TeV: 2 < \rjj\ <5, rjjp rj j2 < 0, p T j > 50GeV, Mjj > 1000GeV. (3.36b) 

We perform numerical fits to the total cross section for —1 ^ r ^ 1 and —1 ^ x ^ 0.5, and derive 
the following 

a{pp hhjj) 


cr(pp - 
cr(pp 


hhjj) sl 
■ hhjj) 


14TeV 


= (l-x) 2 (l-0.86f+4.8x + 0.59f 2 +16x 2 -4.6fx), (3.37a) 


a(pp hhjj) 

S 


= (l-x) 2 (l-0.47f+4.6x + 0.42f 2 +38x 2 -4.1fx). (3.37b) 


100TeV 

We hnd that the cross section of VBF channel is much more sensitive to x than the other two pro¬ 
cesses discussed above. After implementing VBF cuts, the cross section is dominated by longitudinal 
weak boson scattering, and the amplitude has E 2 enhancement which greatly improves the signal 
sensitivity to x in pp(100TeV) collisions. In Fig.[4j we present the cross sections by red curves at 
the LHC (14TeV) and pp(100TeV) collider. These cross sections are normalized to the NLO SM 
prediction D21 at (r, x) = (0, 0). The total cross sections become comparable to that of the gluon 
fusion production over large negative x region, but their dependence on r is weaker. 

In Fig.JTJ we present the distributions for the leading p T of Higgs boson (first row), the di- 
Higgs invariant-mass M ,. (second row) at the LHC (14TeV) (first column) and pp(100TeV) collider 
(second column). In comparison with top-pair associated production of Fig.JhJ more signal events are 
populated in the high p T and M hh regions for x ^ 0 , which is notable even at the LHC (14TeV). To 
further ensure the perturbative expansion of the present effective theory, we will take into account 
the unitarity constraint. We require signal events to obey the conservative bound M hh > A / U1 in 


(2.29). For y/s = 14TeV (lOOTeV) collisions, this allows 84% (31%) signal events under x = —1. 
and 97% (62%) signal events under x = 0.5 . 


10 For the ratio between the VBF signal cross sections in Eq. ( |3.37| ), we note that the QCD A'-factors are largely 
cancelled out and thus this ratio is very insensitive to the A-factors. 
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Figure 7: Parton level distributions of pp —> V*V*jj —> hhjj for the leading p T of Higgs boson 
(first row), the invariant-mass M hh (second row) at LHC (14TeV) (first column), and pp(lOOTeV) 
(second column). In each plot, we set r = 0, and input x = (—1, 0, 0.5) which correspond to 
(blue, red, green) curves. 


4. Full Analysis of gg — > hh — > 6677 at pp(lOOTeV) Collider 

In this section, we study di-Higgs production via gluon fusion by performing a full analysis (including 
Delphes3 fast detector simulations) at the pp(100TeV) collider. We will focus on the gluon fusion 
process gg —>• hh —» 6677 . We construct four kinds of benchmark points, and study the sensitivities 
to probing different regions of the parameter space of cubic Higgs interactions via this channel. Our 
analysis extends the previous Snowmass study |6j by including non-SM-like derivative cubic Higgs 
coupling via model-independent dimension-6 effective operators. We also present a full background 
study which further includes jet-faking-photon backgrounds and contributions of jj 77 due to mis- 
tagging b or b. These improve the analysis of Ref. [6] . 


4.1. Full Simulations for Signals and Backgrounds 

For the present study, we generate the signal and background events by using Madgraph5 and 
Pythia6.2 packages |42] |48|. which are then passed to Delphes3 for detector simulations [49] . 
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Table 1: For signal and background processes, this table presents a X Br, generated events, selected 
events, acceptance, and the expected events at pp(lOOTeV) collider with an integrated luminosity 
of 3 ab _1 . 


Samples 

a x BR (fb) 

Generated Evt 

Selected Evt 

Accept 

Expected 

/i(66)/i( 77 ) (SM) 

3.53 

100000 

3955 

0.040 

418.8 ± 6.6 

666 .( 77 ) 

50.49 

99611 

78 

0.00078 

118.6 ± 13.4 

^( 66 ) 6 ( 77 ) 

0.8756 

68585 

378 

0.0055 

14.5 ±0.7 

tth( 77 ) 

37.26 

63904 

67 

0.0010 

117.2 ± 14.3 

tt'yy 

335.8 

150654 

1 

6.6 xlO -6 

6.75 ±6.7 

tt 'y 

108400 

285787 

0.013 

4.7 xlO -8 

15.2 ±3.2 

6677 

5037 

763962 

11 

1.4 xlO " 5 

217.6 ±65.6 

bbJn 

8960000 

1119406 

0.0051 

4.6 xlO -9 

123.6 ±31.9 

3311 

164200 

813797 

0.056 

6.9 xlO -8 

33.9 ±3.8 

Total background 

- 

- 

- 

- 

647.3 ± 76.0 

s/Vb ( s/vb+s ) 

- 

- 

- 

- 

16.5 (12.8) 


We show the full list of backgrounds in Table]!} All background processes include up to one 
extra parton with MLM matching to avoid double-counting. We do not include bbjj background, 
since after all selection cuts it is negligible compared with other faked backgrounds. The detector 
responses are based on the current performance of ATLAS and CMS. The 6 -tagging operation point 
is chosen to have 75%, 18.8%, and 1% for bottom, charm, and light flavor jets in the central region 
(E t > 50GeV and |r/| < 2.5), respectively. The photon identification efficiency is about 80% for 
photons with E T > 50 GeV and |ry| < 2.5 . For the jet-faking-photon background, we assign a faking 
probability of /■ = 0.0093 exp(—E T /27) as a function of E T (in GeV) of the jet, and scale the jet 
energy by 0.75 ± 0.12 as the photon energy [50]. The mass resolution is 2 GeV for h — > 77 and 
17 GeV for h —> bb at M h = 125 GeV. To be consistent with the signal, we select two tagged 6 -jets 
and two isolated photons in the final states, where each object is required to have E T > 25 GeV and 
\r]\ < 2.5. 

We further impose the mass-window cuts on the invariant-masses of two photons and two 6 -jets. 
Compared with the previous study [ 6 ] , we will narrow down the diphoton invariant-mass window as 
122 GeV < Af 77 < 128 GeV. This would kill another 40% backgrounds beyond the previous case 
with 10 GeV diphoton mass-window. For two 6 -jets, we still impose 85 GeV < M bb < 135 GeV. 

Fig.[ 8 ]shows the normalized distributions of the p T and the sub-leading E T of two selected photons 
(or 6 -jets) in the first two rows. The last plot of Fig. [ 8 ] depicts the reconstructed di-Higgs invariant- 
mass for both signals and backgrounds. Here we only show the representative backgrounds. 

The distributions of faked 6677 and jj 77 are similar to 6677 , while tt 77 and ttrf have too few 
events after selection. For illustration, we present distributions for the SM and two other cases with 
new coupling inputs (r. x) = (—1, 0.5) and (f, x) = (1, —1). We find that including the new 
couplings (r, x) does not significantly change kinematic distributions after full simulation for the 
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gluon fusion production, as we have expected from the parton level analysis in Sec. [3} Hence, for the 
rest of selections, we use the same kinematical cuts as in the Snowmass study [6j. 

We summarize these cuts as follows, 


• Invariant-mass cut: M, r > 300 GeV; 

OO77 

• A R cuts: Ai? 77 < 2.5 , A7?^ < 2.0 ; 

• p T cuts: p T [ 7], p T [b\ > 35 GeV, p T [ 77], p T [bb\ > 100 GeV; 

• Decay angle of h —> 77 in the hh rest frame: | cos 6 h \ < 0.8 J^] 

• Total number n of jets, photons and leptons are required to be n <7 in each event. 

We present the expected signal and background event numbers at ^/~s = 100 TeV and for an 
integrated luminosity C = 3ab _1 in Table[l] For the SM Higgs self-coupling of (r, x) = (0, 0), we 
find the expected signal events to be 418.8 . The expected yield of total background events is 647.3, 
with the largest contributions coming from 6677, bbj'y, bbh( 77) and tth( 77). The resultant signal 
statistic significance is about 16.5 a. With some relaxation of kinematical cuts, we find that the 
sensitivity becomes a bit worse due to increased background contributions, but the overall picture 
remains the same. We have also compared our study with the recent analyses of 6677 channel at 
pp( 100) TeV in the literature [51] [25] . Ref. [51] studied this channel for the SM cubic Higgs coupling, 
and estimated 179 signal events with 447 background events after all cuts and for the same luminosity. 
Our study gives 418.8 signal events and 647.3 background events. The difference is likely due to their 
more conservative assumptions for the detector performance, especially the photon identification 
efficiency, which is lower than ours. In the future, it would be helpful to directly compare the results 
by using the same assumptions for detector performance. Ref. [25] estimated S/\^B = 15.2 under all 
cuts and the same condition, which is in good agreement with ours. 

For the signal analysis, we perform full simulations for parameters within the range —1 ^ r ^ 1 
and — 1 ^ x ^ 0.5. We find that the number of selected signal events can be fitted by similar 


functions as in Eq. (3.34). Under the above cuts, we deduce 


a 


cr Q 


All 


= (l-x) 2 (l -0.55r + 3.4z + 0.11f 2 + 3.9x 2 - 1.2fs) 


(4.38) 


Compared with the parton level fit (3.34b), we see that the cross section becomes less sensitive to 
the parameters (r, x). This is what we would expect from the contamination of parton shower, 
hadronization, and detector simulation. 

To further discriminate r and x dependence, we can utilize distributions in different recon¬ 
structed di-Higgs invariant-mass bins [271 [ 25], which include different kinematic features of contri¬ 
butions from f and x. To efficiently suppress the background, we choose Mhh(= M b j; 77 ) bins as 


11 The decay angle 6 h is defined as the angle between one of the h directions in the di-Higgs rest frame and the 
di-Higgs momentum in the lab frame. 
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Figure 8: Distributions of the sub-leading E T [sub7] and p T [ 77] of selected diphotons for the sig¬ 
nal/background events are presented in the first row. The distributions of Fly [sub 6-jet] and p T [bb] of 
the selected 66 jets are depicted in the second row. The invariant-mass distributions of the selected 
7766 events are plotted in the third row. 
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follows 


M hh bins (GeV): [300,500], [500,700], [700,900], [900,1100], (4.39) 


We note that for the 6677 final state, due to the small branching fraction of h —> 77 and the fast 
decline of gluon parton distribution function, the probe of M hh is not much higher than 1 TeV 
even at the pp(100TeV) collider. Since the derivative cubic Higgs coupling brings in more energy 
enhancement, higher M hh bin is more sensitive to ?. This can be seen from event fits in each bin 
as follows, 


a 


O'sm 


bin 1 


<7 

^sm 

a 

^sm 

<7 

^sm 


bin 2 

bin 3 

bin 4 


(1- ?) 2 ( 1 - 0.82? + 3.4? + 0.17? 2 + 3.3 ? 2 
(1- ?) 2 (1 - 0.42?+ 3.3? + 0.06? 2 + 3.8 ? 2 
(1-?) 2 (1 -0.14?+3.5? +0.04 ? 2 + 5.6 ? 2 
(1-?) 2 (1 -0.03?+4.0?+ 0.03 ? 2 + 8.6 ? 2 


1.5??), 

(4.40a) 

0.95 ??), 

(4.40b) 

0.85??), 

(4.40c) 

0.65 ??). 

(4.40d) 


With increasing M hh , the coefficients of ? terms decrease, while ? terms become more important. 
In passing, we clarify the difference of our analysis from Ref. [251 • The paper [25] simplifies the 
computation by doing hadron-level analysis for the SM case only, and infers the signal rate at other 
points by parton-level analysis with rescaling of hadron-to-parton ratio for the SM, i.e., they assumed 
that the hadron-to-parton cuts efficiency remains the same over the parameter space. We test this 
assumption with our full analysis in the ? — ? parameter space. We find that it works well in 
lower M hh bins, but would induce 0(10%—100%) deviations in high mass binsj^j For the inclusive 
rate, it is not a problem since it is dominated by low mass bins. But, it could affect the conclusion 
of exclusive analysis (cf. Sec. |4.2[ ). For later convenience, we summarize the numbers of selected 
background events for each bin in Tablc[2| 


4.2. Probing Cubic Higgs Interactions via Parameter Space (?, x) 

In this section, we analyze the probe of ? — ? parameter space at the pp (100 TeV) collider with a 
sample data from 3 ab ^ 1 (30 ab _1 ) integrated luminosity. As mentioned in Sec.[ 2 j due to interferences 
with other possible dimension -6 operators, the measurement of single Higgs productions cannot 
uniquely constrain ?. Hence, it is important to independently probe the parameter space of ? via 
di-Higgs production, which receives energy enhancement from the derivative coupling induced by 
2 ■ To study sensitivities to different regions of the (?, ?) parameter space, we choose four kinds 

12 Dcpending on the luminosity, there could be 0(10%) statistical uncertainty in the last Mhh bin at 30 ab -1 . But 
the statistical uncertainties in other bins are much smaller. 
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Table 2: Selected events in different M hh bins for the SM signal and backgrounds at the pp (lOOTeV) 
collider with an integrated luminosity of 3ab _1 . 


M hh bins (GeV) 

[300, 500] 

[500, 700] 

[700, 900] 

[900, 1100] 

h{bb)h( 77 ) (SM) 

200 

170 

52.5 

11.1 

66 / 1 ( 77 ) 

67.1 

31.9 

15.8 

3.81 

Z(66)/i(77) 

11.2 

2.77 

0.46 

0.04 

tth( 77 ) 

97.5 

15.9 

3.22 

0.58 

tt'y'y 

5.41 

1.1 

0.24 

0.0 

tt'y 

13.9 

1.09 

0.16 

0.05 

6677 

188 

23.7 

5.25 

0.32 

66.77 

107 

11.8 

3.44 

1.32 

ini 

30.3 

2.58 

0.82 

0.24 

Total Backgrounds 

521 

90.8 

29.4 

6.37 


of benchmark points, 


Benchmark A : 
Benchmarks Bi, B 2 : 
Benchmarks Ci, C 2 : 
Benchmarks Di, D 2 : 


( ? ) ®)sm (®) 0) 1 

(f, x) = (0, 0.2), (0, 0.5); 

(f, x) = (-0.5, 0), (0.5, 0); 

(f, x) = (-0.5, 0.2), (0.5, -0.5). 


(4.41) 


Benchmark A corresponds to the SM Higgs boson, and the sensitivity in this case can be directly 
translated into a bound on the effective cutoffs of dimension -6 operators (O# 2 , 0<$ 3 ). We use 
Benchmarks Bi and B 2 to represent the cases as predicted by nonminimal coupling model with 


r = 0 and x > 0 (cf. Sec. 2.2). Benchmarks Ci and C 2 correspond to nonzero r and vanishing 
derivative cubic Higgs coupling x . The last two benchmarks Di and D 2 denote the general cases with 
both r and x nonzero. For all non-SM benchmarks, we choose (r, x) values corresponding to the 


effective cutoffs A 2 , A 3 > 500 GeV. Note that the effective cutoff scale A ■ = A/y/$j is not exactly 
the mass scale of an underlying new particle (as the dimensionless coupling ■ could be larger 
than one and is usually less than about 4 - 7 r). One example is the model of Higgs-gravity interactions 
in Sec. 2.2 with dimension -6 operators ( |2.21 )-(2.22). From the viewpoint of effective theory, the 
major issue is to make sure that the energy scale is within the perturbative unitarity bound, so 
that the perturbative analysis is valid. In Fig.[2j the plots (b)-(c) show that for the effective cutoff 
A 2 > 500 GeV, the unitarity bounds on the scattering energy are well above 1 TeV. This justifies our 
perturbative analysis of signal events with di-Higgs invariant-mass M hh <1.1 TeV. 


For each benchmark, we first analyze the sensitivities in different M .. bins as defined in (4.39). 
For a given set of (r, x), the 68 % C.L. contour is defined as follows, 

A Stf, x) 


V B i + S i(r, x) 


= 1 


(4.42) 
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Figure 9: Sensitivity to 5r — 5x plane around Benchmark A, (r, x) = (0, 0). In each plot, the 
dashed (solid) curve depicts 68%C.L. contour with 3ab _1 (30 ab -1 ) integrated luminosity, and the 
dotted line denotes the degenerate direction around the origin. Plots (a)-(d) present the results for 
each Adi 
respectively. 


hh bin. Plots (e) and (f) show the inclusive sensitivity (4.43) and exclusive sensitivity (4.44), 


where signal S i and background B i in Table[2] denote the numbers of selected events in a given bin 
A'lj :^, and ASt(r, x) = 15)(r + 6r, x + dx) — S)(r, x)|. The dependence of signal on the parameters 


(r, x) is determined by the numerical fits in Eq. (4.40). Around the origin of (r, x), it is well ap¬ 
proximated by the linear expansion, S t ~ q + a^r + b^x . It means that the signal is only sensitive 
to the combination afir + b^x , but not the perpendicular direction b t Sr — a t 5x . We call the later 
as “degenerate direction”, along which the signal remains constant nearby the origin. Using the fit 


(4.38), we further derive the sensitivity contour with inclusive data. 

*) 


+ x)] 


= 1. 


(4.43) 


Finally, to fully utilize the information of different M hh bins, we can derive the combined contour at 
68% C.L., 



(4.44) 


which we will call “exclusive” sensitivity, 
which only uses the total rates. 


This is stronger than the “inclusive” sensitivity (4.43) 
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Figure 10: Exclusive sensitivity contours (68%C.L.) in A 2 — A 3 plane for Benchmark A at the 
pp(100TeV) collider with an integrated luminosity of Sab^ 1 (dashed curves) and 30 ab^ 1 (solid 
curves). The two red and blue contours correspond to x 2 x 3 > 0 and x 2 x 3 < 0, respectively. The 
region on the right-hand-side of each contour (and above it) is allowed. 


In Fig. [9] we analyze the sensitivity for Benchmark A, which corresponds to taking the central 
values (f, x) = (0, 0) as in the SM. We present the 68% C.L. contours for each M hh bin in plots (a)- 
(d). Then, we show the inclusive sensitivity contour (4.43) in plot-(e), and the exclusive sensitivity 
contour ( |4.44 ) in plot-(f). For each plot, the dashed (solid) curve depicts 68% C.L. contour with 
3ab _1 (30 ab -1 ) integrated luminosity, while the dotted line shows the degenerate direction around 
the origin. The slope of dotted line varies for different bins of M hh . It is clear that higher M hh bins 
are more sensitive to x, as also noted before [271125]. However, the final sensitivity of a given bin 
also depends on the number of selected events in this bin. Due to suppression in the tail region of 
M,, distribution, event number in the highest bin (purple) could be quite small. This is the case with 
3ab _1 data in Fig.^d), where the sensitivity to x is much lower than that in other bins. Hence, the 
inclusive sensitivity is mainly determined by the first two bins. For 30ab _1 data, there are enough 
events in the last bin to probe x with a good accuracy. Impressively, since various bins are sensitive 


to different combinations of r and x , the exclusive analysis (4.44) makes a big improvement of the 
sensitivity, as shown in Fig.[9)(f ). Note that the exclusive analysis does not improve much of the 
sensitivity for each parameter alone, but helps to break the degenerate direction in the 2-dimensional 
plane. This demonstrates the important role played by the derivative cubic Higgs coupling x in 
the di-Higgs production. It means that gluon fusion production could probe both (f, x) to a good 
accuracy. This is a new point. For comparison, we derive the sensitivity to each parameter alone by 
fixing the other parameter to its SM value. We find that the exclusive sensitivity to Sr is about 13% 
(4.2%), and that to Sx is about 5% (1.6%), for the 3ab _1 (30ab -1 ) integrated luminosity. 

In Fig. 10, we present the exclusive sensitivity contours (68% C.L.) in A 2 — A 3 plane for Bench¬ 
mark A at the pp(100TeV) collider with an integrated luminosity of 3ab _1 (dashed curves) and 
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Figure 11: Sensitivity contours in 5r — Sx plane for Benchmark Bi with (f, x) = (0, 0.2) as shown 
in plot-(a), and for BenchmarkB 2 with (r, x) = (0, 0.5) as shown in plot-(b). In each plot, the 
dashed (solid) curve depicts 68%C.L. contour with 3ab _1 (30 ab -1 ) integrated luminosity, and the 
dotted line denotes the degenerate direction around the origin. The blue and red contours depict the 
inclusive sensitivity (4.43) and exclusive sensitivity (4.44), respectively. 


30ab” 1 (solid curves). The region on the right-hand-side of each contour (and above it) is allowed. 
The cases of x 2 x 3 >0 (if < 0) and x 2 x 3 <0 (if > 0) are shown by the two red and blue con¬ 
tours, respectively. For each contour, the asymptotically flat or vertical behavior gives the sensitivity 
to one operator (when the other is absent), which can be read from the intersection of 68%C.L. 
sensitivity contour in Fig.[9](f ) with each axis. The sensitivities of probing the two operators are 
comparable, A 2 ,A 3 > 1 TeV with 3ab _1 , and A 2 ,A 3 > 2TeV with 30 ab” 1 . For the blue contours, 
the cusps correspond to the end points of ellipse long axis in Fig.[9][f). These cusp regions give the 
weakest 2d sensitivities, A 2 ,A 3 > 0.75 TeV for 3ab 1 data, and A 2 , A 3 > 1.4 TeV for 30 ab 1 data. 
For the red contours, the 2d sensitivity is always stronger. 


In Fig. 11, we present the inclusive sensitivity (4.43) and exclusive sensitivity (4.44) for Bench¬ 
mark Bi [plot-(a)] and BenchmarkB 2 [plot-(b)] by blue and red contours, respectively. The dashed 
(solid) curve depicts 68%C.L. contour with 3ab” x (30 ab -1 ) integrated luminosity, and the dotted 
line denotes the degenerate direction around the origin. The Higgs gravitational interaction predicts 
r = 0 and x > 0. As shown in plots (a) and (b), the sensitivity contours, including slope of the 
degenerate direction, strongly depend on the explicit value of x . Fig.[l2] demonstrates the sensi¬ 
tivities for Benchmark Ci [plot-(a)] and Benchmark C 2 [plot-(b)], where x = 0 and two nonzero r 
values take opposite signs. We find that their shape and sensitivity range are quite similar to that 
of Benchmark A (the SM case). This is expected given the fact that the di-Higgs total cross section 
and invariant-mass ( Mhh ) distribution are much more sensitive to x than r. 


In Fig. 13, we present the inclusive sensitivity (4.43) and exclusive sensitivity (4.44) for Bench¬ 
mark Di [plot-(a)] and BenchmarkD 2 [plot-(b)] to illustrate the features for (r, x) both nonzero. 
Benchmark Di represents the case that the signals in the first two bins of M hh are quite insensi- 
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Figure 12: Sensitivity contours in 5r — 5x plane for Benchmark Ci with (r, x) = (—0.5, 0) as shown 
in plot-(a), and for BenchmarkC 2 with (r, x) = (0.5, 0) as shown in plot-(b). In each plot, the 
dashed (solid) curve depicts 68%C.L. contour with 3ab _1 (30 ab -1 ) integrated luminosity, and the 
dotted line denotes the degenerate direction around the origin. The blue and red contours depict the 
inclusive sensitivity (4.43) and exclusive sensitivity (4.44), respectively. 




Figure 13: Sensitivity contours in 5r — 5x plane for Benchmark Di with (r, x) = (—0.5, 0.2) as 
shown in plot-(a), and for BenchmarkD 2 with (r, x) = (0.5, 0.5) as shown in plot-(b). In each plot, 
the dashed (solid) curve depicts 68%C.L. contour with 3ab -1 (30 ab -1 ) integrated luminosity, and 
the dotted line denotes the degenerate direction around the origin. The blue and red contours depict 
the inclusive sensitivity (4.43) and exclusive sensitivity (4.44), respectively. 
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tive to 5x at the linear order, and the shape of the 68 % sensitivity contour is mainly determined 
by quadratic terms. Although the last two bins still have strong dependence on x, the inclusive 
sensitivity is determined by the first two bins (due to their large rates) with parabola-like shape. 


The exclusive sensitivity is largely improved, especially with 30 ab 1 data. Fig. 13 b) presents the 


sensitivity contours ( 68 %C.L.) for BenchmarkD 2 , where all M hh bins have strong dependence on 
x. The sensitivity to x is significantly enhanced as compared to other Benchmarks P*1 Since the 
sensitivity has little change among different bins, the 68 % contour is only slightly improved by the 


exclusive analysis (4.44). 


In summary, the qualitative feature of sensitivity contours in the 6r — 5x plane can vary sig¬ 
nificantly for different benchmarks. In some cases (such as Benchmarks A, Bi, Ci, C 2 , and Di), 
the exclusive analysis of different M hh bins makes big improvements. In particular, it can break the 
possible degenerate direction around the origin, and impose much stronger constraints on the 2 d 
parameter space even with the di-Higgs production measurement alone. For some other cases (such 
as Benchmarks B 2 and D 2 ), the parameter-dependence of signals appears quite similar in different 
bins. Thus, both the exclusive and inclusive analyses give comparable sensitivities. 

5. Conclusions 


Despite the LHC Higgs discovery, the Higgs boson self-interaction is fully untested so far. It is the 
key ingredient of Higgs potential, and plays vital roles for electroweak symmetry breaking, vacuum 
stability, electroweak phase transition, and Higgs inflation. This is a most likely place to encode new 
physics beyond the standard model (SM). 

In this work, we studied the probe of cubic Higgs interactions via di-Higgs production at hadron 
colliders. We parametrized the new physics of Higgs self-interactions in terms of model-independent 
dimension -6 effective operators in section|2j We take the nonminimal Higgs-gravity interaction as 
an explicit example to motivate such effective operators. The contributions of the two dimension -6 


operators (2.16) to cubic Higgs couplings have different kinematic structures as shown in Eq. (2.10) 


They give different kinematic distributions in various di-Higgs production channels due to the different 
energy-dependence. This is demonstrated in Figs.[5]-[7]of section[3j We also analyzed the weak boson 
scattering and tt scattering at high energies, and derived perturbative unitarity constraints on the 
parameter space in Fig. [2] Among the three channels of di-Higgs production, top-pair associated 
production and vector boson fusion (VBF) production are more sensitive to the energy-enhancement 
in high energy collisions, though their cross sections are generally smaller than the gluon fusion 
production (Fig. [4]). 

In section[4j we performed systematical Monte Carlo analysis of di-Higgs production gg —> hh 
in the decay channel hh —> 6677 by using Delphes 3 fast detector simulations. We computed both 


Note that the plot range of 8x in Fig. 13 b) is much smaller than that in Fig. 13 ’a). 
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signals and full SM backgrounds at the pp(lOOTeV) collider with a Sab^ 1 integrated luminosity, 
as summarized in Tablc[l] and Fig. [8} This channel shows a good potential of discovering cubic 
Higgs couplings in pp(lOOTeV) collisions. Our derived significance is in main agreement with the 
literature j 25j , while a difference from m appears due to the different assumptions about detector 
performance. We further studied the probe of new physics effects in the r — x parameter space with 
full simulations. Since different bins of the di-Higgs invariant-mass M.. exhibit distinctive kine- 
matical features, we used them to discriminate the two dimension-6 operators. We did an exclusive 
analysis to incorporate such kinematical information and obtained a big improvement of sensitivity. 


We further identified four kinds of representative benchmarks (4.41) for the parameter space of cubic 


Higgs coupling, which have qualitatively different features. For each benchmark, we studied the 


sensitivity to the 2d parameter space of (r, x), via both the inclusive analysis (4.43) and exclusive 
analysis (4.44). For comparison, we used two sample integrated luminosities (3ab _1 and 30ab~ x ) of 


the pp(100TeV) collider. For Benchmark A (SM case), the exclusive analysis breaks the degeneracy 
in the 2d plane and makes it possible to probe both (f, x) to a good accuracy by the di-Higgs 
measurement alone. This is demonstrated in Fig.|9] [l0| For one-parameter analysis, we found that 
with a 3ab _1 (30ab _1 ) integrated luminosity, the exclusive sensitivity to r and x are about 13% 
(4.2%) and 5% (1.6%), respectively. Fig. [IT] presented Benchmarks Bi and B 2 with r = 0 and x > 0, 
as motivated by the nonminimal Higgs-gravity interaction. We found that the sensitivity contours 
strongly depend on the size of x . Fig.|12| analyzed Benchmarks Ci and C 2 with x = 0 and different 
values of ? . As expected, the dependence on the change of r is pretty weak. For general regions with 
(r, x) both nonzero, we found that the sensitivity contours behave qualitatively different from the 
SM case of (r, x) = (0, 0), as shown in Fig.|13|(a)-(b) for Benchmarks Di and D 2 . In the case where 
the parameter-dependence of signals in different bins is similar, such as BenchmarkD 2 in Fig.|13[b), 
the improvement of the exclusive analysis over the inclusive analysis becomes rather modest. 


A. Redundancy of Dimension-6 Operators 


In this appendix, we discuss the redundancy of dimension-6 operators in Eqs. (2.4) and (2.5). In the 
SM action, the Higgs sector contains the following terms, 


S sm D 


!dx 4 


(D^H)\D^H) + y 2 H ] H - - y f LHf R + h.c. 


(A.45) 


where L = (f R , /j() T denotes the SU(2)l doublet, and f R the SU(2) L singlet. Then, we can derive 
EOM for the Higgs held, ( D 2 H )t = y 2 H 1 — 2A (H^H)H^ — y^Lf R , and its hermitian conjugate. 
After integration by part, we can rewrite the operator O <^2 as 

20 $j2 = H) = —(WH)d^d^(WH) + (total derivative) 

= [2 (D^H)^(DpH) + H ] D 2 H + (. D 2 H) ] H 

= —20^^ — 2/r^ (H^ H) 2 + 12AO $j3 + ( y^O^j + h.c.), (A.46) 
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where d^(H^H) = H). In the above, we have neglected the total derivative term. We have 

also implemented the SM EOM in the last step, since we only keep operators up to dimension-6. 
With the relation (A.46), we may replace 0$ 4 by other operators, 

1 


A 2 L 


/$, 2 CV 2 + /$, 3^$,3 + 4 + + h.c.) 


M 2 /. 




A 2 


(.WHY + 


/<J> ,2 f< I>,4 


A 2 


O®. 2 + 


/$3 + 6 A / # 4 


A 2 


^$,3 + 


2/<j> j +y ff( t>,4 


2A 2 


+ h.c.^ 


“2 S (^2- a: 4) C, $,2 + Z3 + X4 


3M ‘ to 


4>,3 


+ 


y f. 


Xj + — x 4 ) 0& f + h.c 




(A.47) 


Note that this also shifts the quartic Higgs coupling in the original Higgs potential, but can be 
absorbed by a coupling redefinition, A —y A — /i 2 /<j> 4 /A 2 . At the order of A -2 , the coupling A in 
front of 4 can be replaced by the leading order relation A = M 2 /2v 2 . Hence, for on-shell physical 
amplitudes, we can organize their dependence on (/$ 2) /$ 3 ) /$ 4 : /$ /) via the three combinations 
in Eq. (A.47). 


B. Loop Functions for Triangle and Box Diagrams 

For the analyses of Sec.[3}|4j we need to compute cross sections of the di-Higgs production via gluon 
fusion g{p a )g(p b ) —> h(p c )h(p d ), which invoke loop functions of triangle and box diagrams [52j. The 
triangle loop function is given by 



= Tf [1 + (1 - Tf)f (t f )] , 


(B.48a) 



( -21 
arcsm ,— , 

Vh 

T f > 1, 


f( T f) 

= < 

r ,_ -| 2 

1 1 T f 

j log i-^ " • 

r/ < 1, 

(B.48b) 


where Ty = 4m 2 /s, and s is the partonic center of mass energy. The box loop functions are defined 
as follows, 

[45' + 8m 2 SC ab — 2mjS(S + 2p — 8)(D abc + D bac + D acb j\ 

+m 2 f (2p - 8) [T{C ac + C bd ) + U(C bc + C ad ) - m}(TU - p 2 )D acb ] , (B.49a) 

Gn = s(ru — p 2 ) { m /^ 2 + P~ ~ [SC ab + T(C ac + C bd ) — m 2 STD bac ] 

+ m 2 (U 2 + p 2 — 8U ) [SC ab + U(C bc + C ad ) — m 2 SUD abc ] 

- m){T 2 + U 2 - 2p 2 )(T + U - 8 )C cd 

-2mj(T + U-8)(TU- p 2 )(D abc + D bac + D acb )}, (B.49b) 
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where p = M^/m 2 , S = s/m 2 , T = t/rrtj , U = ii/mj, T = T — p and U = U — p. The two scalar 
integrals are given by 

rd A q 1 


Cij = 


Fijk — 


17r 


(q 2 —m 2 ) ytq+Pi) 2 -m 2 ^ | {q+Pi+Pj) 2 ~m 2 

d A q 1 


17T 


(q 2 —m 2 ) yiq+Pi) 2 -m 2 J ^(g+Pj+Pj) 2— |^(g , +p i +p J +p A .) 2 -m^ 
In the low energy limit s <C m 2 , the loop functions behave as 


(B.50a) 

(B.50b) 


T a = - + 0 [ — 


. m 


2 > 


F u = -- + 0 , G n = C7 


, rrit 


, m. 


In the high energy limit m 2 <C s , they take the asymptotical forms, 


F a = - 


rn. 


1 2 


rrif 

log —r~ + i7r 
s 


m. 


m 4 


m. 


+ 0 -F , F u = 0 -T , G u = 0 -+\ 


(B.51) 


(B.52) 
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